Anharmonic Torsional Stiffness of DNA Revealed under Small External Torques 
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DNA supercoiling plays an important role in a variety of cellular processes. The torsional stress 
related with supercoiling may be also involved in gene regulation through the local structure and 
dynamics of the double helix. To check this possibility steady torsional stress was applied to DNA in 
the course of all-atom molecular dynamics simulations. It is found that small static untwisting signif- 
icantly reduces the torsional persistence length (It) of GC-alternating DNA. For the AT-alternating 
sequence a smaller effect of the opposite sign is observed. As a result, the measured It values are 
similar under zero stress, but diverge with untwisting. The effect is traced to sequence-specific 
asymmetry of local torsional fluctuations, and it should be small in long random DNA due to com- 
pensation. In contrast, the stiffness of special short sequences can vary significantly, which gives 
a simple possibility of gene regulation via probabilities of strong fluctuations. These results have 
important implications for the role of local DNA twisting in complexes with transcription factors. 

PACS numbers: 87.14.gk 87.15.H- 87.15. ap 87.15.ak 



The double helical DNA in living cells is subjected to 
a constitutive unwinding torque created by special en- 
zymes. This forces DNA to fold in a supercoiled state 
similarly to a flexible rod with bending and twisting elas- 
ticity. The supercoiling is long known to play an im- 
portant role in a variety of cellular processes Its 
magnitude changes regularly during the cell cycle and 
in response to environmental conditions, which is ac- 
companied by activation or suppression of certain genes 
0. In E. coli, relaxation of the superhelical stress si- 
multaneously alters activity of 306 genes (7% of the 
genome), with 106 genes activated and other deactivated 
[3J. The genes concerned are functionally diverse, widely 
dispersed throughout the chromosome, and the effect is 
dose-dependent. These and many similar observations 
suggest that the DNA supercoiling is used as a univer- 
sal transcriptional regulator Q, but the corresponding 
physical mechanisms are not clear. 

Detailed studies indicate that the promoter sensitivity 
to supercoiling stems from the recognition of promoter 
elements by RNA polymerase, and that it does not re- 
quire DNA melting or transitions to alternative forms [J] . 
The supercoiling torque is distributed between twisting 
and writhing so that the untwisting of the double helix is 
estimated as 1-2% which is below the thermal noise 
and too small for reliable recognition. However, the ac- 
tion of the torsional stress can be conveyed through a 
property rather than the structure of the double helix. 
The behavior of the supercoiled DNA is governed by the 
interplay between the local bending and twisting fluctu- 
ations. If the bending flexibility or the torsional stiffness 
of the double helix vary with forced untwisting, parame- 
ters of thermal fluctuations could be noticeably affected 
already for short DNA stretches involved in the recogni- 
tion. This idea is appealingand it is supported by some 
earlier data for long DNA [6j-|8(. Local torsional fluctua- 
tions are likely to be involved in regulation directly. In 



bacterial promoters, the optimal linker between the -10 
and -35 elements involves 16 base pair steps (bps), but in 
promoters sensitive to supercoiling it is usually one step 
shorter or longer 0, [t| . One step corresponds to rota- 
tion by 34.5°, which approximately equals the root-mean 
square width of torsional fluctuations for the linker. Very 
strong torsional fluctuations of short DNA stretches are 
necessary for activation of some animal promoters [Toj j . 

Local effects of the torsional stress are difficult to re- 
veal experimentally, but they can be probed by all-atom 
MD simulations. New methods were recently developed 
to apply steady forces and torques to short stretches of 
DNA [llj. In contrast to twisting by periodic bound- 
ary constraints and potential restraints used earlier [l2T - 
Il4j the steady stress emulates local conditions of a short 
fragment in a long supercoiled DNA, which makes pos- 
sible evaluation of elastic parameters under very low 
torsional load corresponding to physiological conditions. 
This method captures linear elastic responses as well as 
the twist-stretch coupling effect under small torques cor- 



responding to physiological degree of supercoiling [ll|. 
Here we present the results of the first computational 
study of the elastic parameters of DNA in such condi- 
tions. 



Dynamics of two tetradecamer DNA with AT- and GC- 
alternating sequences, respectively, were simulated in ex- 
plicit aqueous solution using earlier described protocols 
[ll| . For each duplex, nine 164 ns trajectories of all-atom 
dynamics were computed with fixed torque values in the 
range ± 20 pN-nm, which gives about 3 fj,s of simulations 
in total. Three additional trajectories were computed for 
the GC-alternating fragment for verification. Below we 
consider only evaluation of the torsional stiffness. Other 
methods and protocols are described in Appendix. In 
the harmonic approximation the torsional free energy of 
a DNA fragment of length L subjected to external torque 
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where <E> is the overall winding angle, $ T is its equilibrium 
value, It is the torsional persistence length, and kT is the 
Boltzmann's factor. The equilibrium winding varies with 
the torque as 



(2) 



In the course of MD simulations one measures the proba- 
bility distribution P$ for the winding angle of one helical 
turn which, in the limit of infinite sampling, has a canon- 
ical form 



P$ ~ exp 



2L V 



(3) 



The equilibrium winding is estimated as the time aver- 
age (<&)*, and the torsional persistence length l t is ex- 
tracted from the time variance A|$. The potential of 
mean force (PMF) corresponding to any Gaussian distri- 
bution is quadratic, but if the harmonic approximation 
is truly valid, It must be constant with different r. 

The top panel of Fig. [TJ shows variations of $ r corre- 
sponding to Eq. ([2]). All measurements were taken for 
the central 12 bp stretches, with the two terminal steps 
ignored, which gives about one helical turn. The ampli- 
tude of the forced winding is ± 2 %, i.e. about 0.7° per 
base pair. The straight lines shown have the slopes corre- 
sponding to It obtained under zero torque. In the range 
of torques ± 10 pN-nm the points are compatible with 
a linear elastic response (harmonic elasticity). Beyond 
this range the profile remains roughly linear for the AT- 
alternating sequence, but for the GC-alternating duplex 
evident deviations from harmonicity are found. These 
deviations are reproducible and quite strong. If the It 
value were evaluated by Eq. ([2]) using $ r for r = ±20 
pN-nm it would be about 200 nm. 

The measured torsion persistence length changes with 
the applied torque as shown in the bottom panel. The 
GC-alternating sequence exhibits strong anharmonicity, 
with the twist increase of 1.4° per bps accompanied by 
30% growth in l t . For the AT-alternating sequence, the l t 
profile is nearly flat with a small decreasing trend. This 
trend becomes more visible with stronger twisting (article 
in preparation). The bending stiffness varies somewhat 
beyond the estimated statistical errors, but without reg- 
ular trends. 

Fig. [2] shows the probability distributions P$ for the 
GC-alternating sequence for three representative values 
of r. All of the distributions are close to the analytical 
Gaussians defined by Eq. ([3]) with different l t . Since 
the width of the bells changes, the neat shapes of the 
computed distributions are not due to the harmonicity 
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FIG. 1 Color online. Representative torque dependences ob- 
tained by all-atom MD simulations. The results are shown 
for the overall twisting (top panel), the bending persis- 
tence length (/;,, middle panel), and the torsional persistence 
length (bottom panel) of the AT-alternating (<t>o ~ 363.1°, 
red squares) and GC-alternating (<t>o ?s 381.8°, black circles) 
sequences. The open circles feature the verification tests. 
The straight lines on the top panel correspond to Eq. J2{ 
with It = 124 nm (solid black line) and It = 145 nm (dashed 
red line). The error bars show statistical errors evaluated 
by the method of block averages (see Appendix). In the top 
panel the symbol size corresponds to maximal errors. 



of the torsional potential. These Gaussian shapes result 
from the central limit theorem of the probability theory 
whatever the underlying potential is. As seen in Fig. 
[3l the single-step twist fluctuations at GpC and CpG 
steps produce wide and skewed non-Gaussian distribu- 
tions strongly different from that predicted by Eq. JTJ) 
(see Appendix). With the temperature around 300K, the 
local DNA dynamics goes far beyond the area where the 
harmonic approximation is valid. However, the torsional 
fluctuations of four consecutive bps already give an al- 
most ideal Gaussian. It can be formally described by Eq. 
(P) and ([3]), but the shape of this bell does not correspond 
to the harmonic approximation of the local free energy. 
The Gaussian profile of fluctuations in long DNA is linked 
with the single-step distributions by a linear growth of 
the variance with the chain length. Consequently, not 
just the apex zones of the skewed distributions in Fig. [31 
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FIG. 2 Color online. The normalized probability density 
obtained with different applied torques. Form left to right, 
the MD results are shown for r= -20, 0, and +20 pN-nm 
by green, red, and blue points, respectively. The solid lines 
exhibit analytical distributions Eq. (|3]l corresponding to the 
measured values of It and <E> r . The lower panel displays the 
same data in semi-logarithmic coordinates. 



but their entire shapes contribute. Therefore, the anhar- 
monicity is significant, but hidden. In addition, the twist 
fluctuations at consecutive steps are anticorrelated and 
partially cancel out. 

The asymmetry of the single-step PMFs is the prob- 
able cause of the variable torsional stiffness of the GC- 
alternating fragment. In the first approximation, the l t 
value is proportional to the second derivative of the PMF 
in the energy minimum (see Eq. ((J)). For an asymmetric 
PMF a decrease in l t may be expected when the external 
torque pushes towards the even slope of the energy pro- 
file. In the GC-alternating sequence both single step dis- 
tributions are left-skewed (see Fig. |3]); so the right-hand 
slope of the PMF is steeper than the opposite one, which 
explains the sign of the trend in l t observed in Fig. [T] 
The nearly flat It profile for the AT-alternating fragment 
can be also rationalized because in this case a strong pos- 
itive skewness of TpA steps is partially compensated by a 
negative skewness of ApT steps (see Appendix) . Prelim- 
inary analysis of other sequences reveals that the strong 
negative skewness of the CpG single-step distributions is 
exceptional (see Appendix) . The homopolymer ApA and 
GpG steps are nearly symmetrical whereas the single-step 
distributions for AG- and AC-alternating DNA indicate 
that they would behave similarly to the AT-altcrnating 
fragment. These conclusions should be yet verified in 
more intensive computations, but we expect that for ran- 
dom DNA the macroscopic torsional stiffness should be 
nearly constant because among the steps with skewed 
distributions positive and negative skewness are equally 
represented. In contrast, for short sequence motives an- 
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FIG. 3 Color online. The probability density P$ for GC- 
alternating fragments of one, two and four bps (from top to 
bottom) obtained with r = 0. The solid red lines exhibit the 
analytical distributions Eq. J3J corresponding to the mea- 
sured values of It and <I>o. On the top panel, the distributions 
for GpC and CpG steps are shown in green (left) and blue 
(right), respectively. 



harmonic effects of both signs are possible. They can be 
very significant because biological systems operate with 
much larger torques than we use here. For instance, the 
binding sites of the phage 434 repressor contain a variable 
4 bp spacer that does not interact with the protein and 
supposedly participates in gene control via the sequence- 
dependent elasticity [Hj]. In the complexed state, this 
spacer is always overtwisted by about 30° (l6| . that is 
ten times the amplitude of twisting in Fig. [TJ 

The experimental bending rigidity of free DNA is char- 



acterized by lb ~ 50 nm [17]. The measured l t values 



vary between 36 and 109 nm depending upon the specific 
methods and conditions [18(. Observations of sequence 
effects are rare [3] , and there are a few reports on the in- 
fluence of supercoiling §1. If we assume that MD over- 
estimates the stiffness of DNA uniformly then the con- 
vergent estimate of It is around 90 nm, close to its value 
in single molecule experiments [2^, [2l[ . The bias can be 
due to the neutralizing salt condition in MD or other fac- 
tors (see Appendix). The nearly quantitative agreement 
between MD and experiment is remarkable because none 
of the parameters used in simulations was adjusted to re- 
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produce the DNA elasticity. We hope, therefore, that the 
detailed microscopic picture provided by MD captures 
the qualitative physical trends dictated by the atom-level 
mechanics of the double helix. Our results argue that, 
under normal temperature, the local DNA elasticity is 
strongly anharmonic. Extrapolation from the apparent 
harmonic behavior of macroscopic DNA is not justified 
despite a good agreement with atomistic simulations for 
chain lengths beyond one helical turn 22, [23|. In ad- 
dition, these computational observations shed new light 
upon some earlier controversial issues. 

According to Fig. [TJ with the helical twist slightly 
shifted from the equilibrium value the sequence depen- 
dence of the DNA elasticity can be significantly changed 
and enhanced. The measured torsional stiffnesses are 
similar without applied torque, but diverge with untwist- 
ing. The deformability of DNA is long considered as 
a possible governing factor in the sequence-specific site 
recognition [l5j |. but this mechanism requires strong se- 
quence dependence of elastic parameters com par ed to 
that observed in experiments with free DNA [19(. As 
we see the properties of the relaxed DNA cannot be sim- 
ply transfcred to supercoiled and/or protein bound DNA 
states. Additional studies are necessary to check if the 
elastic properties of the specific binding sites change un- 
der torsional stress. Its mag nitude may be very large in 
some protein complexes |16| . 

Another debated issue concerns the mechanisms of 
gene regulation via DNA supercoiling 2, ..3] . Many of such 
observations are readily rationalized if we assume that 
the sensitive promoters are regulated via the torsional 
stiffness. Even a slight shift in its value has a dramatic 
effect on the probabilities of strong twisting fluctuations. 
Many transcription factors are designed to bind the dou- 
ble helix at two sites separated by a spacer of several base 
pair steps. They can work as sensors of torsional fluctua- 
tions in DNA. A strong twisting fluctuation may be nec- 
essary for binding such factor or for recognition by other 
proteins of a permanently bound torsional sensor. Fig. 
[5] shows that for fluctuations observable during 164 ns, 
physiological modulations of the torsional stress would 
change the corresponding probabilities by several times. 
For less frequent larger fluctuations the effect would be 
much stronger. One can extrapolate the pattern in Fig. 
[2] to events observable in the millisecond time range, and 
this leads to essentially all-or-nothing switching. 

The external torque shifts the distributions in Fig. [2]by 
changing symmetrically the energies of opposite fluctua- 
tions. If the shape of the distributions does not change 
each pair of curves should give a single intercept between 
the corresponding two apexes. However, if the shifting is 
accompanied by widening, one more intercept should ap- 
pear in the range of large twisting opposite to the torque 
direction. For instance, the negative torque shifts the 
distribution in Fig. [5] to the left, but the simultaneous 
widening raises its right wing and, with very large over- 



twisting, the left curve should go above the other two. 
It is seen in Fig. [2] that the vertical difference between 
between the three plots indeed exhibits a reducing trend 
with large This effect is somewhat paradoxical and it 
qualitatively contradicts the behavior of simple models 
where the torsional energy depends upon a single vari- 
able. Our attempts to reproduce it in discrete wormlike 
chains with anharmonic torsion potentials were unsuc- 
cessful. However, such behavior is possible, in principle, 
due to coupling between different degrees of freedom and 
it requires further studies. 

To conclude, it appears that small external torques can 
significantly alter the torsional stiffness of the double he- 
lical DNA. The effect is sequence-dependent, and, under 
variable degree of supercoiling, different stretches of the 
double helix can become locally softer or stiffer. This 
can represent a versatile mechanism of gene regulation 
via the probabilities of strong twisting fluctuations. 



APPENDIX 
Simulation protocols 

Tetradecamer DNA fragments were modeled with AT- 
alternating (d(AT) 7 ) and GC-alternating (d(GC)7) se- 
quences. The starting states for classical MD simula- 
tions were prepared as follows. The solute in the canon- 
ical B-DNA conformation [24[ was immersed in a 6.2- 
nm cubic cell with a high water density of 1.04. The 
box was neutralized by placing Na + ions at random wa- 
ter positions at least 5 A from the solute. The system 
was energy minimized and dynamics were initiated with 
the Maxwell distribution of generalized momenta at low 
temperature. The system was next slowly heated to 293 
K and equilibrated during 1.0 ns. After that the water 
density was adjusted to 0.997 by removing the necessary 
number of water molecules selected randomly at least 5 
A from DNA and ions, and the simulations were contin- 
ued with NVT ensemble conditions. The temperature 
was maintained by the Berendsen algorithm pBj applied 
separately to DNA, water, and ions, with a relaxation 
time of 10 ps. Simulations with external torques started 
from equilibrated states after a few nanoseconds of free 
dynamics. 

The AMBER98 forcefield parameters (H,^ were used 
with the rigid TIP3P water model [HI . The electrostatic 
interactions were treated by the SPME method [2!| , with 
the common values of Ewald parameters, that is 9 A trun- 
cation for the real space sum and /3 « 0.35. To increase 
the time step, MD simulations were carried out by the in- 
ternal coordinate method (ICMD) , j3C" 



31] with the inter- 



nal DNA mobility limited to essential degrees of freedom 
and rotation of water molecules and internal DNA groups 
including only hydrogen atoms slowed down by weight- 
ing of the corresponding inertia tensors. [32|, [33[ The 
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double-helical DNA was modeled with all backbone tor- 
sions, free bond angles in the sugar rings, and rigid bases 
and phosphate groups. The effect of these constraints 
is insignificant, as was previously checked through com- 
parisons with standard Cartesian dynamics [22I 133 ] . The 
time step was 0.01 ps and the DNA structures were saved 
every 5 ps. All trajectories were continued to accumulate 
2 15 points, i.e. to about 164 ns. To verify the results 
for d(GC)7 three additional trajectories were computed 
with torques r= -20, 0, and +20 pN-nm, respectively 
These computations were carried out in parallel on 32 
processors starting from independent equilibrated states. 
In all simulations the B-DNA conformations were well 
conserved without visible slow trends or accumulated de- 
formations. Fig. [IB shows some standard time plots for 
two representative trajectories. It is seen that the overall 
properties of the double helices remain stable and that 
the DNA structures were well-equilibrated before the be- 
ginning of the production runs. 
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FIG. 1 FIG. QJ : The time evolution of the double helical 
structures during MD simulations with applied torque r= 
-20 pN-nm. The results for the GC- and AT-altcrnating se- 
quences are shown in the left and right columns, respectively. 
The helical parameters are averages for the central 12 bp ob- 
tained by program 3DNA The all-atom RMSDs from 
the canonical A-DNA (red) and B-DNA (black) were com- 
puted for the entire fragment length (14 bp). The analysis 
was carried out for 1000 states equally spaced throughout 
the trajectories and all plots were smoothed with a sliding 
window of 82 ps. 



The choice of the fragment length and sequences is 
consistent with the recent computations 11, 35, 3^| and 
it was dictated by the following considerations. The 
length slightly larger than one helical turn is convenient 
for measuring the elastic parameters of DNA [36J . These 
molecules are homopolymers of AT- and GC-units, there- 
fore, they cannot have distinguished asymmetric struc- 
tures like static bends. True homopolymer DNA duplexes 



have special features and, in free MD with the AMBER 
forcefield, these structures deviate from the canonical B- 
DNA stronger than AT- and GC-alternating sequences 
[37| . The terminal base pairs open rather frequently 
during nanosecond time scale MD, which significantly 
perturbs the flanking DNA structure. Because this dy- 
namics cannot be averaged during the accessible dura- 
tion of MD trajectories, we blocked it by applying non- 
perturbing upper distance restraints as explained else- 
where 



11] 



The statistical convergence of MD sampling 
also suffers from rare transitions of backbone torsion an- 
gles to non-standard states considered as forcefield arti- 
facts 38l 39| . In the present simulations such transitions 
occurred mainly in the terminal base pair steps excluded 
from analysis. No a/7 flips were observed in the middle 
d(CG)6 stretches. A few such transitions that occurred 
in the d(TA)g fragments were documented and discussed 
in our previous report [llj . They increased the statisti- 
cal noise in the measured parameters, but did not cause 
overall structural perturbations. Fig. [2^ shows several 
representative distributions of the backbone torsions in- 
volved in non-canonical transitions (38j . The number of 
a/7 flips in d(TA)g was maximum two per trajectory, 
but some of them were reversed. The distribution with 
the maximal population of non-canonical states shown in 
Fig. [2B corresponds to a trajectory with a single stable 
flip. 




FIG. 2 FIG. [2fJ: Distributions 
sions for three representative 
open circles correspond to the 
suits for d(TA) 6 and d(CG) 6 
The solid black line shows the 
with the largest population of 
+4 pN-nm 



of a, J3, and 7 backbone tor- 
trajectories. Plots shown by 
data in Fig. [TJ5, with the re- 
in red and blue, respectively, 
distributions for a trajectory 
non-standard conformers (r= 
, d(TA) 6 ). 



The sampled conformations of the double helix were 
analyzed by program 3DNA [34|, with only 11 central 
base pair steps considered (d(TA) 6 and d(CG) 6 ). Accu- 
rate direct measurement of the length of short double he- 
lices encounters technical difficulties [lH, |4(j. Therefore, 
the corresponding DNA length was computed as 11*0.335 
nm, that is assuming the length of one step correspond- 
ing to experiment. The error in the attributed length can 
be partially responsible for the systematic bias in the ab- 
solute values of other measured parameters, but it does 
not affect qualitative trends. 



External torques 

Steady external torsional load was applied as described 
in detail in our recent report [111 ]. This method dis- 
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tributes forces over selected groups of atoms and com- 
pensates them by reactions applied to other atoms so as 
to zero the total external force and torque. Because the 
forces are applied at different points internal stress and 
deformations are introduced that correspond to overall 
twisting. The method was thoroughly verified in Brow- 
nian dynamics simulations of calibrated discrete worm- 
like chain models [HI]. Notably, it was checked that the 
torque values in the range of interest cause linear elastic 
responses in perfect agreement with theoretical predic- 
tions and negligible side effects. 



Evaluation of statistical errors 

Evaluation of errors in MD simulations is based upon 
the following assertions from the probability theory. Con- 
sider a random variable x with expectation Mi = £ and 



variance Da; 
compute 



and 



a . We can take n samples of x and 



-(xi 



X2 



1 " 
+ X n ) = - V 
n — ' 



X k 



k=l 



s 2 = 



1 ^ 

fc=l 



(x k - x) z 



called the sample average and variance, respectively. 
Both x and S 2 are random variables, with Mi = £ and 



MS* 



a , i.e. x and S provide unbiased estimates of 



£ and a , respectively. It is also known that 

a 2 

Di = — 

n 

and, if x is a Gaussian random variable, 
& _ o_4 



T>S 2 = 2o*/(n-l). 



(4) 



(5) 



Eq. (HJ) and ([5]) are used for evaluation of statistical er- 
rors. 

In our MD simulations the random variable is the wind- 
ing angle of one helical turn, with expectation M$ 
and variance D<F The torsional persistence length is 
computed as 

It = L/D$ 

where L is the DNA length. It can also be obtained 
from the shifts in M$ caused by external torques of dif- 
ferent magnitude. The true values of M<I> and D<& are 
estimated by using, respectively, the sample average and 
variance computed over all n points saved during an MD 
trajectory. However, Eq. (jU and ((5|) cannot be applied 
straightforwardly because they are valid only for statis- 
tically independent samples, i.e. the time intervals be- 
tween the MD states must be suitably large compared 
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FIG. 3 FIG.[3]3: Analysis of statistical errors by the method 
of block averages [4l| . The upper panel displays the results 
for the nine MD trajectories of the GC-alternating fragment. 
For comparison, the lower panel presents a similar analysis of 
eight equivalent BD trajectories of the same DNA fragment. 
The equivalence means that MD and BD trajectories have 
identical durations and saving intervals. The BD simulations 
were carried out under zero torque by using a discrete WLC 
model from earlier reports 11, 36]. The black dashed line on 
the lower panel displays the results for a single eight times 
longer BD trajectory. 



to the torsional relaxation time. The data saving inter- 
val is commonly much smaller, therefore, the errors are 
evaluated by using the method of block averages (4l| . 
The trajectory is successively divided in n' — 2, 4, 2 15 
stretches (blocks) and the sample variances S' 2 are com- 
puted by using n' block averages instead of individual 
samples. When the blocks are longer than the torsional 
relaxation time the block averages are independent and 
S' 2 /n' w const = cr 2 /h, where n is the effective number 
of independent samples provided by the trajectory. This 
value should be used in place of n in Eq. d4j and ((5]). In 
practice, it is convenient to draw the plots of 

n a 

with respect to log 2 n' . When statistical independence is 
reached, such plots exhibit a plateau with J]wn/n = r c , 
which gives the required estimate of h. Parameter r c 
is the effective correlation time measured in trajectory 
saving steps. 

Analysis of the MD data for the GC-alternating DNA 
is shown in Fig. [3J3. For benchmark comparison, we also 
present a similar treatment of Brownian dynamics (BD) 
trajectories of a discrete wormlike chain (WLC) model 



used in our recent studies [llj, |36(. With n' decreasing, 
the plots display emergence of a plateau and the growth 
of statistical noise. From these plots the t c values are 
estimated as 50 and 30 for MD and BD, respectively, in 
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agreement with the corresponding relative rates of tor- 
sional relaxation 36[. The r c w 50 gives n=655 and the 
relative error of 5.5% in the measured l t values. This 
accuracy is sufficient for our purposes. The lower panel 
of Fig. [3J3 also displays the improvement that might be 
obtained with longer trajectories. As expected, with 2 3 
longer BD trajectory the plateau is less noisy, but its 
value is similar. 



Sequence-dependent distributions of twist 
fluctuations 
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FIG. 4 FIG. I4J5: Normalized single-step probability densi- 
ties of torsional fluctuations in MD. The red lines show the 
Gaussian distributions corresponding to the harmonic WLC 
model with the corresponding values of It- The Gaussians 
were shifted to the maxima of the computed distributions. 
The MD data shown are from the trajectories with zero ap- 
plied torques. 



Fig. compares the observed single step distributions 
of torsional fluctuations in d(AT)7 and d(GC)7 with the 
corresponding distributions in equivalent coarse grained 
models, i.e. discrete WLC models with elastic param- 
eters identical to those in MD. The numbers shown in 
the upper right corners display the corresponding mode 
skewness computed as 

(mean — mode) /(standard deviation). 

Similar results for other base pair steps are shown in 
Fig. [5]3. These data were obtained by MD simulations 
of tetradecamer fragments with the corresponding alter- 
nating or homopolymer sequences. Duration of all tra- 
jectories was about 8 ns. 
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